arXiv:1509.06959v2 [cond-mat.mes-hall] 1 Feb 2016 


Perturbation theory for an Anderson quantnm dot asymmetrically attached to two 

snpercondncting leads 

M. Zonda/ V. Pokorny,^’^ V. Janis,^ and T. Novotny^ 

^Department of Condensed Matter Physics, Faculty of Mathematics and Physics, 

Charles University in Prague, Ke Karlovu 5, CZ-12116 Praha 2, Czech Republic 
^Institute of Physics, Academy of Sciences of the Czech Republic, 

Na Slovance 2, CZ-18221 Praha 8, Czech Republic 
® Theoretical Physics III, Center for Electronic Correlations and Magnetism, 

Institute of Physics, University of Augsburg, D-86135 Augsburg, Germany 
(Dated: February 2, 2016) 

Self-consistent perturbation expansion up to the second order in the interaction strength is used 
to study a single-level quantum dot with local Coulomb repulsion attached asymmetrically to two 
generally different superconducting leads. At zero temperature and wide range of other parameters 
the spin-symmetric version of the expansion yields excellent results for the position of the 0 — tt 
impurity quantum phase transition boundary and Josephson current together with the energy of 
Andreev bound states in the 0-phase as confirmed by numerical calculations using the Numerical 
Renormalisation Group method. We analytically prove that the method is charge-conserving as well 
as thermodynamically consistent. Explicit formulas for the position of the 0 — tt phase boundary are 
presented for the Hartree-Fock approximation as well as for its variant called Generalized Atomic 
Limit. It is shown that the Generalized Atomic Limit can be used as a quick estimate for the position 
of the phase boundary at half-filling in a broad range of parameters. We apply our second order 
perturbation method to the interpretation of the existing experimental data on the phase boundary 
with very satisfactory outcome suggesting that the so far employed heavy numerical tools such as 
Numerical Renormalization Group and/or Quantum Monte Garlo are not necessary in a class of 
generic situations and can be safely replaced by a perturbative approach. 
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I. INTRODUCTION 

In the last decade advances in the fabrication of nano¬ 
devices enabled to connect quantum dots with super¬ 
conducting (SC) leads forming superconducting quantum 
dot nanostructures generalizing the conventional Joseph¬ 
son junctions^. Many experimental realizations of this 
concept using various BCS materials for the supercon¬ 
ducting leads {Al, Pb, or Nb) and a great variety of 
quantum dots formed in semiconducting nanowires^*^ 
or dotS)^ carbon nanotubeS)^^— or even single Ceo 
molecules^! demonstrate the versatility of such setups. 
Major advantage of the superconducting quantum dots 
over conventional microscopic Josephson junctions lies in 
the possibility of a detailed control of their microscopic 
parameters, e.g., by tuning the onsite energy by a gate 
voltage. Such a high tunability is promising for potential 
applications of these hybrids in the nanoelectronics (e.g., 
as a superconducting single-electron transistor) or quan¬ 
tum computing as well as for detailed studies of their 
non-trivial physical properties. 

These include Josephson supercurrent, Andreev sub¬ 
gap transport and the way they are influenced by the 
zero-dimensional nature of the superconducting quan¬ 
tum dots with finite-size quantized levels and poten¬ 
tially strong effects of the local Coulomb interaction 
leading to strongly correlated phenomena such as the 
Kondo effect In many cases the system can be very 
well described by a simplest single impurity Ander¬ 


son model (SIAM) coupled to BCS leads^ which, de¬ 
pending on particular parameters, may exhibit so called 
0 — TT transition signaled by the sign reversal of the 
supercurrentThe 0 — tt transition is in¬ 
duced by the underlying impurity quantum phase tran¬ 
sition (QPT) related to the crossing of the lowest many- 
body eigenstates of the system from a spin-singlet ground 
state with positive supercurrent (0-phase) to a spin- 
doublet state with negative supercurrent (7r-phase) i^^— 
In single-particle spectral properties this transition is 
associated with crossing of the Andreev bound states 
(ABS) at the Fermi energy as has also been observed 
experimentally 

A number of theoretical techniques have been used to 
address the 0 — tt transition and related properties of su¬ 
perconducting quantum dots. A very good quantitative 
agreement with the experiment s^^ can be obtained 
in a wide range of parameters using heavy numerics such 
as numerical renormalization group 
and quantum Monte Carlo (QMC)j ^^'^^i^^d^ However, 
both NRG and QMG are time and computational re¬ 
sources consuming. Alternative (semi)analytical meth¬ 
ods based on various, often quite sophisticated, pertur¬ 
bation approaches either around non-interacting limit 
{U = 0) such as the mean-field theory^i^Siiiiii 

slave particleSj ^^d^ and functional renormalization group 
(fRG)^i^ or around the atomic limit (T —^ 0 or A —>• 
0^^ 36,47,48 been used for qualitative and in some 

limits even a quantitative description of the supercon- 
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ducting quantum dot properties. Yet, none of the 
mentioned methods with the exception of the mean- 
field/Hatree-Fock (HF) approximation are sufficiently 
simple and at the same time versatile to serve as a 
generic (semi)analytical solver. HF approximation has 
the attributes of the generic method^ yet, it suffers 
from fundamental conceptual problems, namely, it iden¬ 
tifies the 0 — TT transition with the point of breaking of 
the spin-symmetric solution and attributes the tt phase 
to the magnetic solution of the self-consistent HF equa¬ 
tions. This unphysical breaking of the spin symmetry to¬ 
gether with the ensuing discontinuities of various physical 
quantities even at non-zero temperatures contradicting 
the experimental observations disqualify the unrestricted 
HF approach as a reliable solver for the superconducting 
SIAM. 

Surprisingly, with the exception of a few fragmented 
precursorsit has not been noticed until very 
recently^Si^ that the resummed perturbation theory in¬ 
corporating second-order dynamical corrections to the 
spin-symmetric HF solution yields at zero temperature 
a nearly perfect description of the 0 phase for symmet¬ 
ric leads in a wide range of parameters. The aim of this 
work is to demonstrate that second-order perturbation 
theory is an efficient and reliable method not only for the 
symmetric leads, but also for a more general and realis¬ 
tic case of asymmetric tunnel coupling to different leads 
(i.e., with various values of superconducting gaps). This 
method is numerically much less expensive than NRG 
or QMC. Note that in the general case one has to deal 
with the two-channel Anderson model, therefore, the in¬ 
troduced second-order perturbation theory can be 10 or 
even 100 times faster (depending on parameters and used 
CPU cores) than the fully convergent NRG calculations. 
Simultaneously, this method gives nearly perfect results 
for the physical quantities in the 0-phase at zero tem¬ 
perature in a wide range of parameters corresponding 
to weakly- and intermediately correlated regime, where 
the conventional deployment of NRG is unnecessary. As 
known from previous studies^Si^I the ground-state in this 
regime is the BCS singlet in contrast to the strongly cor¬ 
related regime where the ground state is the Kondo sin¬ 
glet. We illustrate this in Fig. [T] which depicts the ground 
state phase diagram in the U — A plane for quantum dots 
with symmetric leads at half-filling. The crossover region 
between the BCS and Kondo singlets is approximately 
plotted as a gray stripe. The BCS singlet regime covers a 
broad range of parameters (note the logarithmic scale on 
the vertical axis) where second-order perturbation the¬ 
ory is in a nearly perfect agreement with NRG. Thus, we 
advocate this method to be the generic first-choice solver 
for the properties of the 0 phase. 

In order to support this standpoint we carefully ex¬ 
amined formal properties of the approximation such as 
charge conservation, gauge invariance, and thermody¬ 
namic consistency and showed that it preserves all these 
important requirementsThen, we systematically stud¬ 
ied the zero-temperature 0-phase quantities in a wide 



U/F 

FIG. 1. (Color online) Sketch of the phase diagram in the 
U — A plane of the superconducting single-impurity Ander¬ 
son model with symmetric leads at half-filling. Full red line 
separates the singlet and doublet ground states and the shade 
region signals the crossover between the two kinds of singlet 
ground states. The perturbative approach presented in this 
paper works well in the whole BCS singlet regime as demon¬ 
strated in detail in Figs. |6] and [3 below. 


parameter range, paying a special attention to the posi¬ 
tion of the phase boundary between the 0 and tt phases. 
We identified the limits of applicability of our method 
by direct comparison with NRG data obtained via the 
NRG Ljubljana open source codej^i^ At small enough 
temperatures, the boundary depends only weakly on 
temperature^!^ and, therefore, our zero-temperature re¬ 
sults are directly applicable to the existing experimen¬ 
tal data. We finally compared our predictions with two 
experiments with excellent agreement, further justifying 
our claims. 

The outline of the paper is as follows. In the next 
Sec.iniwe introduce the superconducting SIAM, the Mat- 
subara Green function methodology of the perturbation 
theory together with the Josephson current, and ABS for¬ 
mulas in the first subsection, while in the second part we 
study charge conservation, gauge invariance, and thermo¬ 
dynamic consistency conditions to be obeyed by approx¬ 
imations. In Sec. uni we introduce and analyze proper¬ 
ties of the Hartree-Fock approximation and its dynamical 
corrections. In the following Sec. IIYI after a brief sum¬ 
mary of technical issues concerning the evaluation of the 
approximative equations, we present results for the posi¬ 
tion of the 0 — TT boundary first for the case of identical 
leads (equal SG gaps) and then for different leads with 
unequal gaps. Finally, in the last subsection we discuss 
in details the applicability and limitations of our method 
as demonstrated on various single-particle quantities in 
the 0 phase, such as ABS energies and/or induced SG 
gap. In Sec. El we present comparison of quantitative re¬ 
sults of our theory with two existing experiments on the 
position of the 0 — tt boundary. We conclude our work 
in the last Sec. IVII Supporting technical calculations for 
the HF boundary and charge conservation in the second- 
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order perturbation theory are deferred to Appendixes lAl 
andlBl 


II. THEORY AND METHODS 

A. Model and observables 

We use a single-impurity Anderson model^ii^i^iii^^ as 
a model of the quantum dot with well-separated en¬ 
ergy levels connected to two asymmetric superconducting 
leads)^ The Hamiltonian of the system is given by 

H=Hdot+ ^ (Kead+^?)- (la) 

a—R,L 

The first term represents a single impurity with the level 
energy £ and the local Coulomb repulsion U : 

^dot — ^ + Ud\d^dld^, (lb) 

cr=t.l 

where d\ (d^) creates (annihilates) an electron with the 
spin cr on the impurity. The second term of Eq. (Hal) is 
the BCS Hamiltonian of the leads 


^lead — E Aq,E(c _j^^-|-H.C.), 

ktr k 

(Ic) 

with a = L,R denoting the left and right leads and ^ 0 - 
(cakcr) creating (annihilating) conduction electrons. Fi¬ 
nally, the hybridization term between the impurity and 
the contacts is given by 

'Ht = 


All the studied quantities can be expressed with 
the help of the impurity one-electron (imaginary 
time/Matsubara) Green function (GF), which is a 2 x 2 
matrix in the Nambu spinor formalism 


/ 

V 


/ Ga{T - t') , g_-a{T - t')\ 
ySair - t') , G-aij - t’)) 

/ {T[d 4 r)dUr')]) , (TK(r)d_.(r')]) \ 

\(T[fiL^(T)4(T')]) , mdl^{T)d_a{T')])) 
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where the bar denotes the hole function. Since we only 
consider spin-symmetric solutions throughout the whole 
paper, we skip the spin index and we also set e = h = 1 
from now on. The exact form of the unperturbed {U = 0) 
impurity GF can be written as a function of Matsub- 
ara frequencies ojn = (2n -|- l)7r/,5 (see Appendix A of 
Ref. [13): 


Go^iojn) — 
where 


+ s(ia;„)] - £ , A<i,(ia;„) 

AKiWn) , iuJn[l + s{iuJn)\ + Sj 

(3a) 


3(*W„) = E 


ct—L.R 




- 


(3b) 


is the hybridization self-energy due to the coupling of the 
impurity to the superconducting leads. We have denoted 
by Fl^/j = iiPL,R the normal-state tunnel-coupling 
magnitude (with pl,r being the normal-state density of 
states of the respective lead electrons at the Fermi en¬ 
ergy) and 


A$(ia;„) 


E 

a—L^R 


a/A 2 -t UjI 


(3c) 


The impact of the Coulomb repulsion U > 0 on 
the Green function is included in the interaction self¬ 
energy matrix S(*w„) = ( 5 (E).’ 

full propagator in the spin-symmetric situation is deter¬ 
mined by the Dyson equation G~^{iuJn) = G(C^(ia;„) — 
S(za;„). Symmetry relations for the spin-symmetric ver¬ 
sion of the Green function, Eq. ©, reformulated in the 
Matsubara representation^ (Sec. 9.3.3) are G(ia;„) = 
-G*{iuJn) = -G{-iu}n) and g{iujn) = f/*(*w„) = 
g*{—iuJn)- The same applies to the self-energies, i.e., 
E(ia;„) = -T,*{iu}n) = -S(-iw„), 5(ia;„) = S*{iu}n) = 
S*{—iuJn)- Here the asterisk stands for time inversion 
being complex conjugation in the Matsubara formalism. 
Consequently, the interacting Green function explicitly 
reads as 


G{iuJn) = - 


D{iUJn) 


iWn[l + s(iiOn)]+e+ Yi*{iuin), -A$(iw„)-I-5(iw„) \ 

-A|(ia;„) -I- S*{iijJn), *a;„[I -|- s(ia;„)] - £ - E(iw„)y 


( 4 ) 
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The existence and the position of the ABS are determined by zeros of the negative determinant of the inverse Green 
function (i.e., poles of the Green function) 


= - det[G = uj^[l + s{iu}n)f + |£ + T,{iujn)f + |A$(za;„) - S{ioJn)f 

= D{-iuJn) = D*(iuJn) > 0 . 


( 5 ) 


The determinant analytically continued from Matsub- 
ara to the real frequency axis is real within the gap and 
can go through zero D{uJo) = 0 determining the (real) 
in-gap energies ±wo of the ABS symmetrically placed 
around the Fermi energy lying in the middle of the gap. 
ABS are crucial for transport of the Cooper pairs through 
the quantum dot because they usually provide the dom¬ 
inant contribution to the dissipation-less Josephson cur¬ 
rent through the impurity. Furthermore, their crossing 
at the Fermi energy determines the phase boundary be¬ 
tween 0- and TT phases as their zero energy is equivalent 
to the degeneracy of the two lowest-energy many-body 
eigenstates of the system, which is the point of the im¬ 
purity QPT.Si 

The Josephson current out of the dot into the respec¬ 
tive reservoir is defined from the Heisenberg equa¬ 
tion of motion for the particle number in the reservoir 


Ja = 




)M = 






n 


and can be evaluated as a Matsubara sum of the anoma¬ 
lous Green function 




/3 

X Tr 




Oi ' 






^ 2 F„A, 

/3 ^ 


F„A„ 




= -y 

R 


0 

[f/(*w„) 

\Q (^^n) 


G{iuJn) 


- Q(iuJn)e' 




( 6 ) 


where a = L, R as before. For any approximative treat¬ 
ment such as our perturbation expansion in U there al¬ 
ways arises an important question of charge conservation, 
i.e., whether Jl = — Jr and thermodynamic consistency, 
i.e., whether Josephson current calculated by different 
approaches, e.g., directly from Eq. © or as a phase- 
derivative of the associated free energy, gives the same. 
We devote the following subsection to these nontrivial 
fundamental questions. 


B. Charge conservation, thermodynamic 
consistency, and gauge invariance 

Any consistent approximation must respect charge 
conservation, i.e., the Josephson currents through the 


left and right interfaces sum up to zero (due to the 
above convention for the definition of the Josephson cur¬ 
rent as flowing into the respective lead). The condition 
Jl + Jr = 0 can be rewritten with the help of the second 
line of Eq. (|6]) as 


2„ ^ r„A„ 








[Q(iu}n)e - C/*(ia;„)e**“] 


2 

= 'y [/S.%{iLOn)Q{ilOn) “ A$ (iw„ )C/* (lW„)] 
= ^^y[S*{iuJn)GiiuJn) - S{iU}n)G*iiuJn)] 

OJn 

= -^y^S*(iuj„)G(iiOn) , 


( 7 ) 


where we used A^(iciJn) = D(iujn)G(i<^n) + S(iuJn) from 
Eq. (m and reality of D{iujn) from Eq. ([S]). For symmetric 
leads one can choose real A^{iu!n) and consequently the 
self-energy and Green function fulfill S*{iuJn) = S{—ioJn) 
and G*{iuJn) = G{—iuJn)- The charge conservation condi¬ 
tion 0 is then satisfied automatically as can be seen by 
the change of the summation variable —>■ —u>n- Since 
for asymmetric leads one cannot guarantee the reality of 
A$(ia;n) for all frequencies and thus S*{iujn) 7 ^ S{—iuin) 
approximations must be checked for fulfilling charge con¬ 
servation o by explicit verification. 

Apart from a direct approach to charge conserva¬ 
tion via the explicit formula for the Josephson current, 
Eq. ([S]) , we may also employ an indirect one starting with 
the phase-dependent grand potential (“free energy”) of 
the system. The dissipationless Josephson current can 
also be determined as the phase derivative of the ther¬ 
modynamic potential. Approximate calculations of a 
thermodynamic quantity (such as the Josephson current 
here) lead to the same result when different equivalent 
representations are used only in thermodynamically con¬ 
sistent approaches in the Baym sense i^^i^^ 

A thermodynamically consistent approximation can be 
generated from a Luttinger-Ward functional (t>[G]. It 
is represented in terms of the full one-electron Green 
function G, Eq. o, from which the self-energy is de¬ 
termined via a functional derivative. In our case with 
asymmetric leads we have to treat both electron and 
hole variables as independent parameters. Hence, the 
functional derivates determining the self-energies read as 
T,{iuJn) = fd5(j>[G]/6G*{—iuJn) for the normal part and 
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S{iujn) = I3S4>[G]/SG*(—iuJn) for the anomalous one and 
analogously for the self-energies Y,*{iu}n) and S*{iuJn). 
The grand potential then contains both electron and hole 
variables, where the hole variables are decorated with as¬ 


terisks that have the meaning of complex conjugation 
(only) in equilibrium. The grand potential can be repre¬ 
sented with the aid of the Luttinger-Ward functional as 
follows: 


20[G, E] = 0[G] - ^ ^ + G*(-*w„)E(ia;„) + 

-klog[[jw„(l -F s{iujn)) - £ - E(za;„)][ia;„(l -k s{iu}n)) + e + E*(za;„)] - [A<i,(za;„) - 5(za;„)][A(J,(ia;„) - 5*(za;„)]]} . 

( 8 ) 


Complex variables G{iuin),G*(iiOn), E(za;n), S*(ia;„) 
as well as C/(zw„), C/*(za;n), S{iuJn),S*(iujn) are varia¬ 
tional parameters the physical values of which are de¬ 
termined from the stationarity of the grand potential 
r2[G, E]. Due to the electron-hole symmetry the elec¬ 
tron and hole contributions to the thermodynamic po¬ 
tential are identical, hence, the factor 2 on the left-hand 
side. An approximation is thermodynamically consistent 
and conserving if we are able to determine explicitly the 
Luttinger-Ward functional (/>[G]. Such approximations 
are called 0 derivable. 


Reliable and physically acceptable approximations 
should not only be thermodynamically consistent but 
must be gauge invariant. Observables in Josephson junc¬ 
tion setups depend on the phase difference between the 
leads but they cannot depend on the absolute values of 
the two phases. In other words, the physics must be 
invariant with respect to a global phase shift ^l,r t 
A$, which is a manifestation of gauge invariance. 
Obviously, the building elements of the theory, the Green 
functions, Eqs. (153 and (|H), are not invariant and we 
must always check that they enter the measurable quan¬ 
tities, such as the supercurrent (jS]), in a way that pre¬ 
serves gauge invariance. The resulting self-energy E(za;„) 
and consequently G(za;„) transform equally as Go(zu;„) 
(I3al) under the gauge transformation in thermodynami¬ 
cally consistent approximations. That is, only the off- 
diagonal elements pick up the global phase shift (with 
the respective sign) from which we can immediately see 
that the Josephson current ([5]) is indeed gauge-invariant 
as it should be. 


We can also use the functional of the grand poten¬ 
tial r2[G, E] to prove gauge invariance of a (/)-derivable 
approximation. In gauge-invariant theories, thermody¬ 
namic potentials depend only on the phase difference, 
i.e., D[G,E]($i,$ji) = D[G,E]($i - $^z). Conse- 

,, , .■ T - 

quently, charge conservation - = 

■I’r) ^ immediately follows. 


III. APPROXIMATION SCHEMES 


Since the exact expression for the self-energy of the 
superconducting Anderson impurity model is unknown, 
we have applied the standard Matsubara resummed per¬ 
turbation theory in the interaction strength U to sum¬ 
ming up one-particle-irreducible diagrams for the self¬ 
energy in terms of the dressed one-particle Green func¬ 
tion. To avoid the unphysical spin polarization of the im¬ 
purity, which can be easily obtained within the resummed 
(dressed) approach from the self-consistent solution of a 
nonlinear equation for the single-particle Green function, 
we restrict the solution to the spin-symmetric case. This 
results in the situation when for certain parameters at 
zero temperature there exists no solution for the Green 
function. The breakdown of the solution coincides with 
the crossing of ABS energies at the Fermi energy, i.e., 
with the 0 — TT quantum phase transition. Thus, while 
the 0-phase, which can be smoothly connected with the 
noninteracting case U = 0, can be captured by the per¬ 
turbative approach, the tt phase with doubly degenerate 
ground state is beyond the reach of this simple perturba¬ 
tion theory. We explicitly demonstrate this concept for 
the Hartree-Fock solution, where all quantities at QPT 
can be addressed analytically, but the same scheme car¬ 
ries over to higher orders of the perturbation theory, in 
particular to second order, being the main focus of our 
study. The very possibility of a description of the zr-phase 
within a (suitably modified) perturbative approach re¬ 
mains an open question as we discuss in more detail later 
in Sec. II V Gl 








and 
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A. Spin-symmetric (restricted) Hartree-Fock 
approximation 

Mathematical expressions for first-order perturbation 
expansion in U, Hartree-Fock (HF) contributions, read 




where the HF Green function reads 


G 


HF/ 


o = - 


1 f -|- £ -f 






HF" 


— A^(iuJn) -b 5^^ ^ 

*a;„[l -b s{iujn)\ - e - 


, with the determinant (9c) 


D^^{iujn) = [1 -b s{iuJn)f + [e + -b |A$(*w„) - 5^^ 


(9d) 


The hole (asterisks) functions are obtained from com¬ 
plex conjugation of the equations for the electron func¬ 
tions. Obviously, the frequency-independent self-energy 
of the HF approximation neglects any dynamical corre¬ 
lations caused by particle interaction. Nevertheless, it is 
still capable to describe qualitatively the 0 — tt quantum 
phase transition even without the necessity of the com¬ 
mon, yet questionablebreaking of spin-reflection 
symmetry i^°i^^ Therefore, it is a useful demonstration 
tool of the basic ideas concerning the model as well as 
a worthy etalon of more elaborate methods. The explicit 
formula for the HF phase boundary for the completely 
symmetric case = An, F^ = Fj^ was derived in our 
previous work^; here we provide a general solution. 

Before that, though, we discuss charge-conservation 
and gauge-invariance properties of the HF approxima¬ 
tion. If we insert Eq. (I9b|, into the last line of Eq. © we 
see that it is satisfied. Thus, the HF self-energy yields 
a charge-conserving approximation. Similarly, the HF 
self-energy transforms under the gauge transformation 
identically to the (dressed) Green function as it should. 
Furthermore, HF self-energy can be derived from a man¬ 
ifestly gauge-invariant Luttinger-Ward functional 



Gonsequently, HF approximation is both charge conserv¬ 
ing as well as thermodynamically consistent. 

Now, we turn to the calculation of the HF phase 
boundary. The self-consistent HF equations read as 


U ^ - A^iiuJr^) 

/3 ^ D^^iiLUn) ■ 


Further manipulations of the above equations stated in 
Appendix |A] lead to the following set of equations at 
zero temperature (since we are interested in the phase 
boundary) for the auxiliary quantities Ed = -b £ and 

A = V r 


Ed =e+Z-U 



du) Ed 
2tt D^^{iuj) ’ 


OL 



dio 




r,i^a 


Ao 




27r 





( 11 ) 


Glose to the QPT the inverse denominator l/D^^{ioj) 
is dominated by its zeros at the ABS energies ±Wo (i.e., 
£)HF(j_^o) = 0) which become zero at the QPT. There¬ 
fore, we may use the expansion of the determinant to 
the lowest (second) order in uj reading^ D^^{iuj) ~ 
E^ + |(5p — [1 + Ta/^a]^ which gives us for the 

ABS energies close to the QPT wq ~ -b |dp/(l -b 













7 


^^ra/Aa). Obviously, the position of QPT coincides 
with the situation where Ej, = S = 0. Close to the tran¬ 
sition the integrals are strongly dominated by the poles 
at iwo and we can approximately evaluate the first two 

I 


leading contributions in inverse ABS energy (first of the 
order 1 /ojq while the second of the order 1; all other terms 
are at least of order wq and, thus, irrelevant at the tran¬ 
sition) as follows: 


f 

f 


du! 


duj 1 1 

2tT D^'^iiu}) 2{l + J2a^a/^a)'^ J-oo 2wo(l + Ea ra/Aa)^ ’ 


(12a) 


dw f{u 


2tt D^^{iuj) 


f ^ 

Jo ^ 


duj /(w^) 




(12b) 


uJo—Ed—S—0 


for a smooth function f{x) vanishing at zero, i.e., f{x —> 0) = 0. Using these approximations, we finally arrive at 


Ed 


1-b 


1-b 


U 


2/BjTFF(l+EaWAa) 

u 


= s + 


u 

2' 


2/BjTFF(l + EaWAa) 

with B representing the band contribution 




UB, 


B = 


OO 

/ 


dui 




1 - 


Ac 


V'AfW 


0 


1 + Ec 




Ea 


i^Q 




- 1 


In Appendix lAl we discuss the formula (I14|) in more detail in the symmetric case Al = Ah = A. 
To obtain the phase boundary, we sum up squares of the two equations m- 


{E^d + 1^1') 

, U 

2 o 

= Ev) + 

V rae*‘^“ -f UB 

L 2^7/2+ |<5P(l + ^„r„/AEj 

V 2 ; 

a—L,R 


which yields at the phase boundary -|- |i5p = 0 an implicit equation for the borderline 

2 


C/2 


4(l + EcWAa)2 


= e + 


U 


E r.e 

a—L,R 


i^a 


UB 


(13) 


(14) 


(15) 


(16) 


Omitting the band contribution B from Eq. (HU) one gets a very simple approximation of the boundary called gener¬ 
alized atomic limit (GAL)^: 


ETTeEV^" +ri + r« + 2rLr«cos($L-$«). (i7) 


Interestingly, as we show later on, at half-filling {e = 
—U/2) the simple GAL boundary lies typically very close 
to the results obtained via the NRG method and/or the 
second-order perturbation expansion. This is not sur¬ 
prising, since at half-filling the HE approximation repro¬ 
duces the atomic limit exactly. Hence, one can expect 
that GAL, being a generalization of the atomic limit to 
non-integer occupation of the dot, delivers quite reliable 


results near the charge-symmetric state. 

Eurthermore, we may use Eq. (USD for Hnding 
\/E‘^ -j- |5|2 close to the phase boundary. Using the 
implicit-function theorem we see that a solution with 
> 0 exists on one side of the boundary (more¬ 
over, the side containing the noninteracting U —^ 0 limit, 
i.e., corresponding to the 0 phase) while \/E^ |<5P < 0 

on the other side. There, we must conclude, no solution 











































to the restricted HF equations da exists, only when 
one allows for breaking of the spin symmetry (i.e., finite 
magnetization), which is however unphysical for a zero¬ 
dimensional impurity system, the appropriately extended 
HF equations (ITOl) do have a solution^iiii Since we do 
not want to resort to an unphysical symmetry break¬ 
ing to obtain the tt phase, we must conclude that the 
perturbative spin-symmetric solution breaks down at the 
phase boundary as expected from a general conceptual 
viewpoint)^ Although these findings have been explic¬ 
itly demonstrated on the level of the HF approximation, 
they are actually fully general and apply to any order 
of perturbation theory, in particular also to the second 


order which we are going to address now. 

B. Second order: Dynamical corrections 

It was shown in previous studie a^°’^^ that inclusion of 
dynamical corrections beyond the static HF into the self¬ 
energy can dramatically improve the quantitative pre¬ 
dictions for both the position of the phase boundary and 
the physical quantities in the 0 phase. Already, first cor¬ 
rections from second order of the perturbation expansion 
were sufficient to reproduce fairly well the results of NRG 
in the case of identical leads. Second-order contributions 
to the self-energy read as 
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(18b) 


(18c) 


is the two-particle bubble consisting of the normal and anomalous parts and Vm = 27rm//3 is the m-th bosonic 
Matsubara frequency. Analogously to the HF case before we can explicitly verify the charge conservation condition 
©; calculations are more tedious this time and we present them in Appendix [BJ 

The Luttinger-Ward functional of this second-order correction to the Hartree-Fock approximation reads as 


[G{iuJk)G* {-iu}n)G{iuJn + ivm)G*{-iu}k - iVm) + 2G{iu}k)G*{-iu;n) 
xQi^iujji ( iUJk ^^m) T Qi^^kpQ ( ^a^n)f^(^^n “f ( '^^k ^^m)] 



(18d) 


It is manifestly gauge-invariant. These first two orders the one-particle level. The higher contributions to the 
of the perturbation expansion are well controllable on 
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self-energy become more complex and their classification 
demands to introduce two-particle vertices as discussed 
in detail in Ref. [5§. Therefore, we resort just to the sec¬ 
ond order of the perturbation theory, which proves to be 
fully sufficient in the BCS-singlet regime for weak and 
intermediate coupling. The second order self-energy cor¬ 
rections (together with the first-order HF counterparts) 
are inserted into Eq. to obtain a self-consistent non¬ 
linear functional equation for the Green function as a 
function of frequency. Unlike the HF case, the resulting 
equations for the Green function components defy an¬ 
alytical treatment and must be solved numerically. In 
the following we refer to this approach as the full self- 
consistent dynamical correction (FDG) approximation. 

As discussed previously for the symmetric leads^ 
nearly identical results can be obtained in the weak cou¬ 
pling regime by evaluating the dynamical self-energies 
(fTHl) using just the fully converged self-consistent HF so¬ 
lution as the input into the Green function (DC approx¬ 
imation). The DC approach can be represented by the 
following algorithm: 

1. Compute the HF Green function as described in 

Sec, mra 

2. Compute the second-order contributions to self en¬ 
ergy S(^)(ia;„) and 5^^^(ia;„) using formulas (fTSl) 
with G^^{iu}n) and 0^^{iujn) instead of G(ia;„) 
and 0{iu!n). These second-order contributions stay 
fixed throughout further calculations. 

3. Compute the DC self-energies E^®(za;„) = -|- 

T,^'^\iuJn) and 5°^(iw„) = S^^l{iuJn) with 

and 5*^^^ = 5^^ in the first iteration. 

4. Compute the DC Green function G^^{iuJn) using 
definitions o and m with Ti{iujn) = E^'^(iu;„) 
and S{iu!n) = 5^^(*a;„). 

5. Compute the first order contributions to self ener¬ 
gies: E(i) = f and 5(1) = 

6 . Repeat steps 3 to 5 until the convergence criterion 
of the self-consistency is achieved. 

The algorithm implies that the convolutions in the 
second-order self-energies are evaluated just at the begin¬ 
ning of the procedure. The fixed dynamical self-energies 
are then used to calculate self-consistently the first-order 
contributions to the self-energies. 

Note that the DC approximation is numerically more 
stable close to the phase transition than the FDC ap¬ 
proximation (see Fig. [5]). In addition, it allows us to 
study the intermediate coupling regimes, where full self- 
consistent approaches often fail to give physically cor¬ 
rect solution. It is fairly well known that the fully self- 
consistent second-order approximation as well as its ex¬ 
tensions via sums of ladders and chains fail for inter¬ 
mediate coupling of impurity models. They not only 


smear the Hubbard satellite bands^, they also miss the 
Kondo physics^i That is why simplified self-consistencies 
often provide better approximations than the fully self- 
consistent onesi^ 


IV. RESULTS AND DISCUSSION 

We provide a comparison of the ground-state, i.e. 
zero-temperature, results obtained via the perturbative 
method discussed above with those obtained using the 
NRG approach which is a reliable nonperturbative nu¬ 
merical method for the ground state properties;^ For 
NRG calculations we used the NRG Ljubljana open 
source code^i^ mostly with the logarithmic discretiza¬ 
tion parameter A set to the value of 4 as is common for 
double-channel problems. We also tested and used other 
values of A (see Fig. 0^) and found out that for most of 
the studied cases, the phase boundary is not sensitive to 
A. This is in compliance with the findings discussed in 
Ref. m. 

Various theoretical studies showed that the 0 — tt 
phase transition, where the supercurrent changes its sign, 
is accompanied by a smooth crossing of the Andreev 
bound states2i;22ii2iii^ii^ which is in agreement with the 
experimentsAlthough the perturbation approach 
without spin-symmetry breaking can not be easily ex¬ 
tended into the tt phase and, therefore, does not show the 
actual crossing of the ABS;^ the ABS smoothly reach the 
Fermi energy at the border of the 0-phasej ^°i^^ In Fig. 
we plot examples of the ABS dependencies on the phase 
difference $ for symmetric coupling F/, = F/{ (Fig.[2^) as 
well as on the right coupling F n for asymmetric coupling 
Fi = 1.44r/j (Fig.[^). We identify the point where both 
(bound and anti-bonding) ABS reach the Fermi energy 
with the boundary of the 0-phase, i.e., with the point of 
the quantum phase transition. This is fully supported by 
the NRG data. It can be seen in Fig. [2] that the cross¬ 
ing of the ABS obtained by the NRG (bullets) coincides 
with the merger of the ABS obtained via the DC approx¬ 
imation (solid red lines). The dashed blue lines in Fig.[^ 
were obtained via the FDC approximation. One can see 
that the FDC and DC results practically coincide in both 
presented cases apart from the very close neighborhood 
of the phase transition, where the FDC becomes numer¬ 
ically unstable. 

It should be stressed that analytic continuation of the 
Matsubara formalism to the real frequencies is necessary 
for the study of ABS. This usually leads to quite compli¬ 
cated formulas for the Green functionshowever, some¬ 
times this approach is numerically more stable than the 
Matsubara formalism. Nevertheless, the continuation to 
the real axis can be avoided if one is interested only in 
the phase boundaries. The determinant Dlfujn) yields 
zero at the phase transition point exactly for itOn = 0. 
Therefore, one can use the smooth dependency of D(0) 
on different model parameters for the direct estimation 
of the phase boundary. 
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FIG. 2. (Color online) Andreev bound states energies depen¬ 
dence on the phase difference <& for the symmetric coupling 
r_L = Fij (a) and on the coupling F r for asymmetric coupling 
ri = 1.44r_B at <E> = 0 (b). The solid lines were calculated via 
the DC approximation, the dashed lines via the FDC approx¬ 
imation and the bullets were obtained using NRG with A = 4. 
The point where bound and anti-bonding states merge/cross 
at the Fermi level (zero energy) is identified with the 0 — tt 
quantum phase transition. 


A. Phase diagrams for An = Ah 

Setups where An = An = A in practice mean that 
both leads are made of the same material. Since this 
is a common situation in the experiment, this case has 
been intensively studied^ We start our discussion with 
three phase diagrams known from the literature.— The 
main reason for this is to show that the simple second- 
order expansion theory is sufficient for a broad range of 
parameters and even in the regimes where usually much 
more elaborated techniques are used. Simultaneously, we 
compare the results of the three approximations, namely, 
HF, DC, and GAL and discuss their limitations. 

In Fig. [3^, we plot the ground-state phase diagrams in 
the Tn — U plane for <5 = 0 and <I> = tt at half-filling 
(e = —17/2). The left-right lead asymmetry is fixed by 
the coupling ratio F^ = 1.44 Fh. The first thing that 
should be noticed is that the HF approximation without 





e/A 

FIG. 3. (Color online) Phase diagrams in the Tr — U plane at 
the half-filling (a) and for e = —2.5A (b) and in the F—e plane 
for U = 5A (c). We compare the phase boundaries calcu¬ 
lated by NRG with the spin-symmetric HF, the second-order 
PT/dynamical corrections (DC), and generalized atomic limit 
approximation (GAL). 


broken spin-symmetry yields a qualitatively correct de¬ 
scription of the phase boundary at half-filling. However, 
its phase-boundary curve is close to the NRG border only 
for small Goulomb interaction {U < 2A). A much better 
result is obtained by using its simplification, namely, the 
GAL approximation, which neglects the band contribu¬ 
tion. This implies that the HF approximation overesti¬ 
mates the contribution from the bands. Gonsidering its 
simplicity, the GAL method provides a very good and 
fast approximation of the phase boundary at half-filling 
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even for U ^ A. Nevertheless, regarding the accuracy it 
is overperformed by the DC approximation. One can see 
that the DC border reproduces the NRG data (for A = 4) 
almost perfectly. This shows that the perturbation the¬ 
ory with the simplest dynamical corrections can lead to 
a correct estimation of the quantum phase boundary for 
the studied system. This statement is true not only for 
the symmetric-leads case studied in Ref. but also for 
the experimentally more relevant setups. 

A similar agreement between the DC and NRG phase 
boundaries can be seen in Fig. [3 }d, where the phase di¬ 
agrams are plotted away from half-filling at e = —2.5A. 
On the other hand, the bias of the level energy e strongly 
influences both the HF and GAL curves. The GAL bor¬ 
der approaches the NRG data only around U = 5A, i.e., 
just near the half-filling occupation. The HF border is 
way off in the whole plotted range. In addition, both HF 
and GAL results drift away from the DC and NRG bor¬ 
der with the increasing U. Because of the structure of 
Eq. (ITSl) . the GAL and HF borders approach each other 
for {7 —>-0, where the term UB vanishes, as well as for 
I e -I-17/21 ^ A where the hrst term dominates the right 
hand side of Eq. (HU). The latter one can be seen in 
Fig. [5J:. Here, we plot the phase diagram in the F^; — e 
plane for moderately large Coulomb interaction U = 5A 
and two values of 4). The HF and GAL borders coincide 
far away from half-filling for both values of $. As before, 
despite its simplicity, both these approximations yield 
a fair qualitative agreement with the NRG data. How¬ 
ever, with the exception of the GAL curve near the half¬ 
filling, both approximations fail to reproduce the NRG 
data quantitatively. In contrast, the DC border matches 
the NRG in a broad region around s = —U/2 and even 
outside this region the difference between NRG data and 
DC curve is much smaller than A. This shows that the 
proper treatment of the frequency dependence of the cor¬ 
relation effects is crucial for the quantitative description 
of the 0 — TT transition. 

The most evident (qualitative) discrepancy of the DC 
curve from the NRG are the “humps” at the bottom of 
the phase diagrams Figs. and|3^. This is not surpris¬ 
ing as here U ^ T and therefore we are on the edge of the 
usability of the perturbation expansion in U (see Fig.[T])- 
On the other hand these humps resemble a formation of 
the island-like phase diagram known from the previous 
NRG study by Oguri et al^. They showed that in case 
of strongly asymmetric gaps ^ An and decreasing 
U a re-entrant doublet region appears as an island in the 
Fij — e plane (see Fig. 12 in Ref. IdOll . However, this is 
not the case in Fig. 0^. For the fixed ratio F^/F/j and 
An = An, the humps are only a consequence of the per¬ 
turbative treatment. Unlike in Ref. IdOl they are obviously 
disappearing with the decreasing U (Fig. 0^.) and no is¬ 
land structure appears even for U = A/2. Nevertheless, 
the situation is different if one varies the ratio Tn/Tn- 
In Fig. 0)3 we show that the condition An ^ An used 
by Oguri et al. is not necessary for obtaining the island 
structures of the 7r-phase. All three approximations show 




(e+0.5U)/A 

FIG. 4. (Color online) Phase diagrams in the Fij —e plane for 
different values of U, ^ — n and asymmetric leads Fz, = 2ri{ 
(fixed ratio) (a) and Fn = A (varying ratio Tr/Tl which 
induces the vr-phase island structure) (b). 


these structures when the ratio Vn/^R is varied. Never¬ 
theless, as before, only the DC approximation matches 
quantitatively with the NRG outside the half-filling for 
plotted values of U. 

Because we have observed this behavior also for other 
values of F^/A (for another example see the reentrant 
behavior at An/An = 1 for ?7 = A^ lines in Fig. [5^), 
we conclude that the island structures are primarily a 
function of the difference between the hybridizations F^ 
and Ffl;. In general, the tt phase is destabilized by the 
increasing differences between the couplings. As this can 
be in principle tuned in the experimenti one can expect 
that the island structures can be verified experimentally. 
Regarding the U dependence, one can see that the in¬ 
creasing Coulomb interaction inflates the 7r-phase area 
in both panels of Fig. 01 

We use Fig. 0^ to demonstrate two technicalities. 
First, it is an illustrative way to show that all three ap¬ 
proximations converge to each other with decreasing U 
as it should be. All borders practically coincide for the 
lowest U = A/2. Second, we present a test of the NRG 
for the largest value of Coulomb interaction U = 6A cor- 
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responding to the most correlated case. There is only a 
very weak dependence of the NRG phase boundary on 
the parameter A which makes the generic calculations at 
A = 4 highly reliable. 


B. Phase diagrams for Al ^ Ar 

Significantly less attention has been paid, both experi¬ 
mentally and theoretically, to the general case Al ^ Ar 
than to the identical leads Ar = Ar. This can change 
in the near future as a recent experiment on a carbon 
nanotube quantum dot coupled to Nb fork and A1 tun¬ 
nel probe^^ not only showed that such a setup is tech¬ 
nically possible, but in addition it revealed a nontrivial 
formation of ABS. Theoretical efforts so far have been, 
however, mainly focused on the Ar = Ar case despite 
the fact that previous theoretical studies dealing with 
some special cases of non-identical leads showed that the 
lead difference can strongly affect the quantum phase 

transition . 

We demonstrate this for a wide range of ratios Ar/Ar 
in Fig. Oi. Here, the phase boundary obtained via the 
DC approximation and NRG is plotted in the Tr/Tr - 
Ar/Ar plane for various values of U at half-filling. We 
have set the coupling to the left lead to T^ = Ar/2 and 
the phase difference to $ = tt. One can see that the DC 
curves are in excellent agreement with the NRG results 
in the plotted range of Ar/Ar which spans two orders of 
magnitude (note the logarithmic horizontal scale). This 
once again proves the reliability of the DC approximation 
also for unequal gaps in the leads. 

We can observe two different kinds of phase diagrams 
in Fig. [S^. Two phase boundaries separating the 0- and 
TT-phases are present for 17 = Ar (solid curves), but only 
a single phase boundary is realized in the plotted region 
for U > 2Ar (dashed curves). This is in compliance 
with the study of Ar Ar case in Ref. as well as 
with the opening of island structures as a consequence of 
increasing U observed in the phase diagrams plotted in 
Figs. |1 }d and[S]D. 

Oguri et showed that in the limit Ar —>■ oo the 
model (Hap can be mapped onto a single-channel model 
where the right lead is replaced by an onsite supercon¬ 
ducting gap at the impurity, i.e., the standard super¬ 
conducting atomic limit^^i^l performed for the right lead 
only. Our observation that the critical F/j depends only 
weakly on the ratio Ar/Ar for the large right gap and 
weak Coulomb interaction (see U = Ar boundaries in 
Fig. [S^) justifies applicability of this simplified model. 

The GAL and HF phase boundaries are in a fair qual¬ 
itative agreement with DC and NRG even tor Ar ^ Ar 
case. Nevertheless, both these simple approximations fail 
quantitatively in positioning the critical curves in the 
phase diagram. This can be seen in Fig. [5j) where the 
0 -phase boundaries are plotted in the Tr/Tr - e plane 
for Ar = Ar/A and F/j = Fj;,/2. On the other hand, 
the DC approximation largely coincides with NRG even 




FIG. 5. (Color online) a) Critical tunnel-coupling ratio 
A'r/Tl as a function of SC-gap ratio Ar/ Ar for Fi = Al/2, 
<I> = TT at half-filling, b) Phase diagrams in the Tr/Vr — e 
plane with asymmetric lead SC gaps Ar = Ar/ A calculated 
via NRG and DC approximation for U/ Ar = 1, 2, 4 and com¬ 
pared for U = 4Al with HF and GAL approximations. 


for U = AAr and fully reproduces the 7r-phase island 
structure obtained with NRG for 17 = A^. 


C. Applicability and limitations of the method 

Despite its reliability in a wide range of the input 
parameters the present method has natural application 
limits. We can roughly divide them into two classes. 
The first one is related to the very conceptual founda¬ 
tion of the method in the (resummed) perturbation the¬ 
ory which implies its breakdown at the quantum phase 
transition and impossibility to reach the tt phase with 
the doubly-degenerate ground state. Any quantities in 
the TT phase are presently inaccessible by our approach. 
The impossibility to take into account the tt phase has a 
serious consequence in that the method cannot address 
nonzero temperatures. At nonzero temperatures both 
the 0 and tt phases coexist and this coexistence results 
in a temperature-smoothened behavior of all quantities 
around the transition. This is, however, completely ne- 
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glected in the present perturbative treatment that is built 
upon only the singlet equilibrium state, 0 phase, even 
if it becomes metastable (at non-zero temperature) and 
should actually be replaced by the spin doublet, tt phase. 
Mixing of the singlet and doublet states is relevant only 
in a close vicinity of the phase transition, since far away 
from it the physics is still described well by either of the 
states. Yet, this conceptual limitation is serious, espe¬ 
cially, since there is no clear way how to circumvent it. 
There is, however, a perturbative expansion for elemen¬ 
tary excitations and Green functions in systems with a 
degenerate equilibrium statei^ 

The second class of limitations concerns the standard 
fact that perturbation methods have limited range of 
applicability given by the level of sophistication of the 
included perturbation contributions. Such theories typi¬ 
cally do not explicitly break down but they become quan¬ 
titatively highly imprecise and eventually useless. Our 
approach is conceptually meaningful in the 0 phase and 
can be used not only for determination of the position of 
0 —TT phase boundaries, but also for calculation of various 
(single-particle) quantities such as the Josephson current, 
QD occupation, proximity-induced local gap, energies of 
ABS etcj^i^ Since it is a perturbation expansion in the 
Coulomb interaction truncated at the lowest-order dia¬ 
grams, it is clear that it ceases to be reliable for large 
enough U. This can happen in various ways depend¬ 
ing on the values of the other model parameters A, T, 
and £. We have already encountered such a situation in 
Figs. [5J: and |3^ where there was an obvious discrepancy 
(although not too severe) for small T close to charge- 
degeneracy points. Small T effectively increases the im¬ 
portance of the Coulomb interaction via the increased 
ratio C/r, pushing the system close to the atomic limit, 
where the complementary perturbation expansion in T is 
a more suitable choicei ^^i^^ 

For decreasing F and fixed A, the perturbation 
expansion around the atomic limit works fine, but 
if one allows also the SC gap A to decrease com¬ 
parably to the (normal state) Kondo temperature 

Tk ^ \J ^ exp (l + ^)], the system enters into 

a strongly correlated Kondo state where any simple per¬ 
turbation theory inevitably fails. We demonstrate this 
crossover from the conventional BCS singlet to the Kondo 
singiet^l and gradual failure of the DC approximation in 
Figs. |6] and 0 for symmetric leads at half-filling. In Fig. [6] 
we present the phase diagram in the A — {7 plane studied 
previously in the literature using different methods in¬ 
cluding the NRG— and the expansion around the atomic 
limit (for A much larger than the characteristic ener¬ 
gies of the dot) based on the self-consistent description 
of the Andreev bound states (SC ABS)— We compare 
the phase boundaries obtained via these methods with 
the DC and GAL boundaries in Fig. |6l Note, that the 
DC approximation is so good that we had to use the loga¬ 
rithmic scale for the A-axis to visualize deviations of the 
DC boundary from the NRG data. The DC boundary 



U/T 

FIG. 6. (Color online) Phase diagram in the A — U plane 
for (p = 0. NRG solution (black bullets connected by dashed 
lines) is compared with the DC approximation (red full line), 
GAL (blue dotted line), and the SC ABS method taken graph¬ 
ically from Fig. 3 in Ref. [3^ (green dot-dashed line). Notice 
the logarithmic scale on the vertical A axis, where the ar¬ 
rows point out the values of the gap used for curves plotted 
in Fig. [3 The grey stripe marks the region where DC results 
start to depart from the NRG in Fig. [7] i.e., the position of the 
crossover from the BCS to the Kondo singlet ground state. 

departs from the NRG points for A/ttF < 0.03, never¬ 
theless, even in this parameter range the DC boundary is 
still much closer to the NRG boundary than the SC ABS 
curve, which can be attributed to the violation of the 
large-A assumption inherent to the SC ABS approach. 
However, the GAL branches off from NRG points signif¬ 
icantly with decreasing A. 

The horizontal arrows in Fig. [5] correspond to the 
values of the gap for which the one-particle quantities 
are plotted in Fig. [7^ (ABS energies wq) and Fig. [Tio 
(proximity-induced local gap A^ = U (d^di); the curve 
A = O.OdF is not displayed just for clarity). There is al¬ 
most a perfect agreement between DC and NRG curves 
for A/F > 0.1 in the entire 0-phase. For smaller values 
of A and sufficiently large U both ABS energies and A^ 
obtained using the DC approximation depart from the 
NRG points with increasing ratio U/T and can even be¬ 
come numerically unstable (dashed lines). We estimated 
the region where the DC curves start to deviate signif¬ 
icantly from the NRG by the grey stripe in Fig. [S] We 
consider this to be the crossover region between the BCS 
and Kondo singlet ground states and, therefore, also the 
edge of applicability of the DC approximation. As we 
will show in the next section, this edge still leaves plenty 
of space for the second-order expansion in U to be the 
method of choice for the description of real systems. 

In Appendix [B] we have proven that the FDC approx¬ 
imation is charge conserving in the general case and that 
DC approximation is charge conserving for the exper¬ 
imentally generic case of equivalent gaps. The charge 
conservation for the DC unfortunately does not extend 
to the general case, giving it the same status as the con¬ 
ventional implementation of fRG— In Fig. [51 we compare 
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FIG. 7. (Color online) Limitations of the DC approximation 
— comparison of the NRG (bullets) and DG approximation 
(solid lines) results for symmetric setups at half-filling. Inter¬ 
action strength U dependence of the ABS energies (a) and 
scaled proximity-induced gap Ad = U {d^di) (b) for the four 
values of the gap A shown by arrows in Fig. |6] and taken 
from Ref. (the curve for A = O.OIF was omitted in the 
lower panel just for its readability). The breakdown of the 
DC method can be seen for large interaction U/V and small 
SC gap A as indicated in Fig. [S] The dashed part of the 
A = O.OIGF lines is numerically unstable. 


the numerical results obtained by these two methods with 
the NRG. The supercurrents at the left/right junctions 
are plotted as functions ofe for U = A^, $ = 7r/2, 
Fi = 2TR = Ai,/2 and three values of Ar/Ar. We used 
the FDC approximation to calculate the supercurrent for 
e+U/2 < 0 and the DC approximation for e+U /2 > 0. It 
can be seen that both approaches are in excellent agree¬ 
ment with the NRG. However, only the FDC approx¬ 
imation fully conserves the current in the general case 
Ar Ar. This is shown in the insets of Fig. [5] where 
the details of the current are plotted close to half-filling. 
The DC approximation conserves current for Ar = Ar 
but not for Ar ^ Ar as illustrated in the right inset. 
The difference between Jr and — Jr is nevertheless very 
small, below 1% for the used parameters (see the vertical 
J-axis scale in the inset), which justifies the widespread 
usage of the DC approximation in this work as a faster 
and numerically more stable (especially around the phase 



(e+U/2)/AL 


FIG. 8. (Color online) Supercurrent through the left/right 
junctions as a function of £ for U = Ar, 4> = 7r/2 , F_l = 
2r_R = A_l/ 2, and Ar/Ar = 1,1/2,1/4. The current was 
calculated using the FDG approximation for e + U/2 < 0 and 
DC approximation for e + U/2 > 0. The insets show details 
close to the half-filling. 

boundary) alternative to the FDC approximation, which 
still gives trustworthy results even for Al yf Ar. 

Finally, we mention yet another aspect of the DC 
method, which is its ability to calculate also spectral 
functions. This is possible by making use of analytic 
continuation of the Matsubara formalism to the real fre¬ 
quencies as shown in Refs. [13, HH, and|13- Here we just 
plot the typical normal and anomalous spectral densi¬ 
ties for the asymmetric coupling to the leads calculated 
using the DC approximation (solid red line), FDC ap¬ 
proximation (dashed-dotted blue line) and NRG (dashed 
black line) in Fig. |9l Discretization parameter A = 4 and 
the logarithmic-exponential broadening of the data with 
the broadening parameter b = 0.15 together with the z- 
trick (see the manual to Ref. Issi) was used for the NRG 
plot. Considering the simplicity of the (F)DC approxi¬ 
mations and the broadening of NRG curves, which is fully 
responsible for the discrepancies around the band-edge, 
the spectral functions are in a very good agreement, no¬ 
tice especially the perfect agreement both in the position 
and weight of the ABS. 

V. COMPARISON WITH EXPERIMENTS 

A. Grenoble experiment [Phys. Rev. X 2, 011009 

( 2012 )] 

Realization of a fully tunable superconducting carbon- 
nanotube quantum dot SQUID by Maurand et al^ al¬ 
lowed to determine the phase diagram for the 0 — tt phase 
transition in the F — e plane experimentally (Fig. 1101) . 
In the experiment the Coulomb interaction was mea¬ 
sured from the finite-bias spectroscopy of the Coulomb- 
blockade diamonds and estimated to be [7 = 0.80 ± 0.05 
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FIG. 9. (Color online) Normal — AG/tt (a) and anoma¬ 
lous — SJCz/tt (b) spectral densities for U = 4A, Fh = A, 
Tl = 2TR (asymmetric coupling), $ = 7r/2, and e = —U 12 
(half-filling) calculated using the DC approximation (red solid 
line), FDC approximation (blue dashed-dotted line) and NRG 
(black dashed line). Discrete Andreev bound states within the 
gap are represented by arrows whose height is determined by 
their weight. 


meV. The superconducting gap A k, 0.08 meV was deter¬ 
mined from the peaks in the nonequilibrium (finite-bias) 
differential conductance. The analysis of the maximum of 
normal-state conductance showed that the tunneling am¬ 
plitudes to the leads were balanced — therefore, the sym¬ 
metric setup (r/2 = ri=rfl;) was assumed. Hybridiza¬ 
tion r for different tunings of the setup was estimated 
from the half-width at half-maximum of the Kondo reso¬ 
nance in the finite-bias conductance. The authors argue 
that the Kondo screening plays a key role for the 0 — tt 
phase transition in their device. This statement is sup¬ 
ported by a quantitative comparison of the position of 
phase boundary with their theoretical predictions based 
on the SC ABS approximation)^^ which is a perturbative 
method based on the superconducting atomic limi t^^’^^ 
with expansion for finite gap A. The key role of the 
Kondo screening is emphasized also from the identified 
operating regime of the experiment marked in Fig. 7 of 



(e+U/2)/U 


FIG. 10. (Golor online) Phase boundary between 0- and tt- 
phase as a function of e and F/[/. Both the experimental data 
(with the estimated Coulomb repulsion U = 0.8 meV and SC 
gap A = 0.08 meV) and the theoretical curve calculated using 
the self-consistent description of Andreev bound states (SC 
ABS)^ were taken graphically from Fig. 6 of Ref. [^. The 
DC phase boundary was calculated for D/A = 10 and <E> = 0. 

Ref. [lO. Considering this and the limitations of the DC 
approach discussed above it is surprising that the actual 
DC phase boundary calculated for the same parameters 
as SC ABS ([7/A = 10) is very close to the experimental 
one as shown in Fig. [TUI Although the SC ABS phase 
boundary is closer to experimental points than the DC 
boundary, that is still within the experimental error bars. 
Therefore, we conclude that the DC method performs 
well even beyond its expected validity range and can be 
applied to a very broad range of real superconducting 
quantum dots. 


B. Orsay experiment [Phys. Rev. B 91, 241401(R) 
(2015)] 

The most recent experimental study2^ of the supercon¬ 
ducting carbon-nanotube quantum dot confirmed theo¬ 
retical predictions that the 0 — tt phase transition can 
be controlled not only by the gate voltage but also by 
the superconducting phase difference $ tuned by the 
magnetic flux piercing the SQUID loop with the carbon- 
nanotube Josephson junction. The superconducting gap 
of the leads A = 0.17 meV and the Coulomb interaction 
U = 3.2 meV, both with uncertainty ~ 10%, were experi¬ 
mentally determined by standard methods (see the previ¬ 
ous subsection). Total hybridization F = F/{-|-rL = 0.44 
meV and the asymmetry Tr/Tl = 4 of the couplings 
were obtained by comparing the measured normal-state 
finite-temperature linear conductance with the one ob¬ 
tained from CT-INT quantum Monte Carlc4^ calcula¬ 
tions for the Anderson impurity model analogously to 
Ref. [ 2 ^ The experimental results were compared to 
QMC calculations and an excellent agreement was ob- 
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e(meV) 

FIG. 11. (Color online) Comparison of phase boundary in the 
<0 — £ plane obtained experimentally with different theoreti¬ 
cal predictions. Both the experimental data and QMC points 
were taken graphically from Fig. 4b in Ref. [2^. Experimen¬ 
tally determined parameters (for details see the main text) 
A = 0.17 meV, U = 3.2 meV, F = Fij -I- Fi, = 0.44 meV, 
and Th/Vl =4 are subject to roughly 10% uncertainty. Ac¬ 
cording to Ref. the QMC curve was horizontally shifted 
to match the experiment. We plot the DC approximation for 
U = 3.2meV (solid red curve), the same result shifted by Se 
to overlap the experiment (dashed red curve), and result for 
U = 3.44 meV within the experimental uncertainty without 
any further modification (blue full line). 


served both for the current-phase relation as well as for 
the shape and width of the 0 — tt boundary in the $ — e 
plane. However, a shift of the energy level Se = 0.28 meV 
of unknown origin was needed to overlap the experimen¬ 
tal data with the theoretical curve. 

In Fig. [TT] we have calculated the 0-phase boundary 
with the DC approximation (solid red line) and compared 
it with the experimental and QMC data taken graphi¬ 
cally from Fig. 4b in Ref. [2^ One can see that after 
introducing a small shift Se = 0.14 meV, exactly as it 
was done for the QMC results, the simple second order 
perturbation theory reproduces (the shape and width of) 
the phase boundary almost perfectly (dashed-dotted red 
curve in Fig. HD). We have taken the advantages of the 
DC approximation (its simplicity and speed) to check 
how the boundary depends on the variance of used pa¬ 
rameters. We have found out, that although the shape 
and width of the boundary are quite robust within the 
10 % uncertainty, the ^-position of the boundary is very 
sensitive to the value of U. No shift of the phase bound¬ 
ary is needed if the value U = 3.44 meV within the U- 
uncertainty range is used instead as it is shown in Fig. 1111 
(blue solid curve). This leads us to the conclusion that 
the deviations between theory and experiment observed 
in Ref. [2^ are within the experimental uncertainty. More¬ 
over, this is yet another demonstration of the usefulness 
of the DC approach for real systems. 


VI. CONCLUSIONS 


To summarize, we presented a detailed study of the 
self-consistent second-order perturbation expansion in 
the interaction strength of the superconducting single¬ 
impurity Anderson model. Based on a thorough analysis 
of its properties, we showed that it can reliably substi¬ 
tute time and resources consuming numerical methods 
such us the NRG or QMC for the study of the 0 — tt 
phase transitions and properties of the 0-phase in super¬ 
conducting quantum dots for a broad range of param¬ 
eters. It can be the method of first choice for realistic 
setups with asymmetric tunnel couplings and even for 
unconventional setups with different SC leads. We dis¬ 
closed its big potential by successful fits of two existing 
experimental data sets for the 0 — tt phase boundary, in¬ 
cluding the suggestion for a plausible explanation of the 
existing discrepancy between the newest experiment and 
corresponding QMC results. 

The approach can be straightforwardly applied to any 
single-particle quantity in the 0 phase such as super¬ 
current, local occupation and proximity-induced super¬ 
conducting gap, or energies and weights of the Andreev 
bound states including the position of the 0 — tt quan¬ 
tum phase boundary at zero temperature. Due to its 
perturbation-theory roots it, however, conceptually fails 
in the description of the tt phase and, consequently, also 
in the description of the finite-temperature properties 
close to the phase boundary. A possible remedy of the 
perturbation approach to reach also the tt phase with the 
doublet ground state remains an open challenge, which 
in view of the successes of the method in the 0-phase is 
worth taking up in future studies. 
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Appendix A: Hartree-Fock phase boundary 


Initial version of HF equations (fTUl) can be recast into 
the simpler form (ED by introducing auxiliary quantities 
Ed = +e and 6 = J2a=L r rae*'^“ —5^^ and further 
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using the general identity 

_ 1 ^ioj„0+ Wn[l + s{iLUn)] _ 1 . . , ^ 

/3^ D{lUn) 2 ^ ^ 

valid for any spin-symmetric GF (j4]) due to a sum rule 
reflecting the fundamental anticommutation relation 


— ^ ' g®‘^riO G{iu)n)\ 


i^„o+ 2za;n[l -I- ^(wn)] -I- T.*{iWn) - S(ia;n) 
l 3 D{iuin) 

G{t -t' ^ 0 ") -h G(t - t' ^ 0 “) =< dU> + < dS >= 1 . 


(A2) 


Large-frequency behavior of the self-energy is 

limited by a constant (in case of the HF approximation; 
otherwise it generically decays as l/iw^), which allows us 
to drop the phase-convergence factor in the above sum 


and using the symmetry relations = S(—iwn) 

(implying real and D{—iu}n) = D{iwn) ([SI) we get 

T.* {iu]n)-'s:{iu>n) _ ^ 

D{iuJn) 2^0Jn D{iuJri) 

thus proving the required identity m- 
Determinant D^^{iuJn) explicitly reads 


{iuJn) = . 


1 + E 

OL 

1+E 


\/A2 -h a;2 


+ [£ + EHF]% 


TaAo, ^ ^jjp 




2 -|- (4;2 

a ' n 


Ef, 


E^o 




( A^ 

\ \/A 2 -I- a ;2 


-1 +(5 


El + + 



(A3) 


Because close to the QPT both E 4 and 6 are close to zero, 
the second term in the brackets multiplying can be 
safely neglected in the calculation of the phase boundary 
(as is done in the main text) since this term is effectively 
of higher order in the a;„-expansion. 

Finally, we discuss the band contribution term B (fTil) . 
Its name derives from the fact that when the integral over 
the Matsubara frequencies (m is Wick-rotated to the 


real frequencies, it only contains the continuous (band) 
part of the spectrum, i.e., it does not encompass any ABS 
contributions. The general formula can be recast into a 
more compact form for the generic case with equal SC 
gaps Ai = /S.R = A. Using the substitution w = Asinht 
and mutually canceling the common 2 sinh^ | terms in 
the numerator and denominator of the integrand, we ar¬ 
rive at the expression 


B = 






^2-71 

Jo 


dt cosh^ t 

("cosht -I- ^ cosh^ ^ + 




. , o 


(A4) 


which generalizes the symmetric case F/, = F/j = r/2 Appendix B: Charge conservation for the dynamical 
(and = — 4*_r = 4’/2) studied previously in Ref. [sO- corrections 

1. FDC 

As a special case of Eq. © we now consider the charge 
conservation in the second-order approximation for which 
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the “current defect” reads 


and the bubble contribution (118c 




/3 

with the anomalous self-energy (I18bp 
1/2 


(Bl) xii^m) = -bT. [Giiiyjn + iuJk)G{iu}k) + G{wm + iuJk)Q*{iijJk)] ■ 

(B3) 


Separating the quantity = 6Jn^ + SJa^^ into two 

S ^{iLiJn) = —^ (B2) parts Corresponding to the normal and anomalous Green 

functions constituents of the bubble, respectively, we get 


ATJ 2 

X! + i^k)G{iu}k)GiitOn + Wm)G*{iUJn), 

^ UJn^tOki^m 

4^2 

X! + i^k)G* {iUJk)G{iuJn + iym)G*{iUJn)- 


(B4) 


<jJri 

Using the symmetry relation G*(iuj) = G{—iuj) we can manipulate the first formula as 
2^/2 

= ^ E [G{Wm + iUJk)G{iuJk)G{iuJn + iVm)G*{iiOn) - G*{iVm + iUJk)G* {iuJk)G* {iuJn + il'm)G{iUJn)] 


OJn ,CJfc jl'rr 


/33 
2^/2 

2^2 

= ^ [G{Wm + iUJk)G{iuJk)G{iuJn + iVm)G*{iiOn) - G{il'm + iU!k)G{iuJk)G*{iuJn - Wm)G{iUJn)] 


= 0 , 


(B5) 


where we have used substitutions ujk —>■ —ujk and Vm —>■ 
—Vm in the second term of the sums between the second 
and the third lines and then the shift of the summation 
variable uin — i^m ^ <^n in the last step. Analogously, the 


anomalous contribution can be simplified with help of the 
substitution Vm —^ —Vm and shift of variables u}n,k 
^n,k + Vm as follows: 



= 0 , 


[G{iVm + iuJk)G*{iuJk)G{iUJn + iVm)G*{iuJn) 


[G{iVm + iuJk)G*{iU}k)G{iUJn + iVm)G*{iuJn) 


G*{iym + iUJk)G{iUJk)G*{iOJn + WTn)G{iUJn)] 


G*{-iVm + iUJk)G{iUJk)G*{iU}n - Wm)G{iU}n)] 


(B6) 


which finalizes the required proof of the conserving na- 2. DC 

ture = 0 of the FDC approximation. 

As we have shown numerically in Sec. II V Cl the DC 
approximation is charge conserving for identical gaps 
Ak = Ak = A. This can be proven analytically by show¬ 
ing that both G^'"{ioJn) and are real which is 

a sufficient condition for = 0 in Eq. m- By making 
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use of the gauge invariance it is possible to introduce a The second-order contribution to the anomalous DC 
global phase shift self-energy is also real because it reads as 


(psh = arctan 


( rL-TR 
{Tl + Tr 



(B7) 




/3 


(BIO) 


such that (pL = (psh — (pR = psh + </'/2 for which 
AApuJn) is real for all frequencies. Consequently, the 
Q^^{iuin) and 5^^ are real too, because the equality 


235 ^^ = 5 "^^ - 5 "^^" = -255 


-HF 


HF* 


/3 ^ 


D^^pUJn) 


(B8) 

can be generally fulfilled only when 55^^ = 0. Assuming 
the contrary, i.e., the existence of the special solution 


-p D^^iiLUn) 

CJji 


(B9) 


implies from Eq. (fTUl) the identity ^ \ ~ 

which would in turn mean that Eq. (l9b| is fulfilled for 
any 5^^. Correspondingly, the same must be true also 
for Eq. (IB9I) which can be easily contradicted, e.g., by 
taking limit 5^^ —>■ oo. Reality of then follows 

directly from the Eq. (1^ . 


where both Q^^{iujn) and the bubble contribution 
(see Eq. (IlScI) ! are real. However, the first or¬ 
der contribution to the anomalous DC self energy reads 
as = E g^'^iiujn) where 

g^'^pujn) = +'5^^^(ia;„) - A^{iu}n)j 

(Bll) 

withall5l^'(ia;ji), D^^{iujn), A^{iujn) being real. There¬ 
fore, using the same argument as for 5^^ and g^^{iujn), 
one can show that 5^^^ is real. Consequently, also 
5°^ = 5^^) -I- 5^^^ and t/°®(ia;„) are real which was to 
be proven. 

Note that this proof does not carry over to the general 
case Ar ^ Ar due to the lack of existence of a global, 
i.e., frequency-independent phase shift to make A^iiuin) 
real for all a;„’s. The DC approximation is thus not con¬ 
serving for non-identical leads as revealed in our numer¬ 
ical results, although the observed current conservation 
breaking is very weak (see Fig. [5]). 
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